function F = solve_eqm_cf1_outer(param, const, eqm)
%% Counterfactual 1: Baseline
 
    % Initial n(q): initial equilibrium
    q_indicator = eqm.q_indicator;
    nq_new = cal_nq(const, q_indicator); 

    % Initial integrals: initial equilibrium
    integ = eqm.integ;

    % Initial Pi(q,E,I)
    eqm_cf.Pi00 = ones(1,param.Q);
    eqm_cf.Pi01 = ones(1,param.Q);
    eqm_cf.Pi10 = ones(1,param.Q);
    eqm_cf.Pi11 = ones(1,param.Q);

    % Settings
    tol = 1e-12;
    maxIter = 100;
    maxHistory = 5;
    maxRepeat = maxHistory+1;
    
    % Placeholder/initialize (check repeat)
    q_opt_index_history = zeros(param.Nf,maxHistory); 
    switch_id_history = zeros(param.Nf,maxHistory+1);
    switch_qindex_history = zeros(param.Nf,maxHistory+1);     
    fix_flag = 0;
    repeat = zeros(1,maxHistory);
    
    % Figure
    figure; 
    
    % Iteration
    err = 1;
    iter = 0;
    while (err >= tol && iter < maxIter)

        % save previous H quality choice
        if iter > 1
            q_opt_index_history = [q_opt_index, q_opt_index_history(:,[1:end-1])];
        end   
        
        % save previous H switch id/qindex
        if max(repeat(2:end)) > 0           
            temp_switch_id = zeros(param.Nf,1);
            temp_switch_id(switch_id) = switch_id;
            switch_id_history = [temp_switch_id, switch_id_history(:,[1:end-1])];            
            temp_switch_qindex = zeros(param.Nf,1);
            temp_switch_qindex(switch_id) = switch_qindex;
            switch_qindex_history = [temp_switch_qindex, switch_qindex_history(:,[1:end-1])];
        end
        
        % new guess
        Pi00 = eqm_cf.Pi00;
        Pi01 = eqm_cf.Pi01;
        Pi10 = eqm_cf.Pi10;
        Pi11 = eqm_cf.Pi11;
        nq = nq_new;
        iter = iter + 1;
                
        % solve equilibrium
        eqm_cf = solve_eqm_cf1_inner(param, const, integ, eqm); 
        
        % choose new quality (grid search)
        [q_opt_index, ~] = grid_search(param, eqm_cf.EPiQ0, eqm_cf.EPiQ1, eqm_cf.ExportProbQ, eqm_cf.ExportCutoffQ);        

        % fix quality after maxRepeat
        if max(repeat(2:end)) >= maxRepeat 
            if fix_flag == 0
                % find id and qindex
                fix_id = nonzeros(max(switch_id_history,[],2));
                fix_qindex_history = switch_qindex_history(fix_id,:);
                fix_qindex = max(fix_qindex_history,[],2); 
                % set flag
                fix_flag = 1;
                repeat = repeat*0;
            elseif fix_flag > 0
                % find id and qindex
                fix_id2 = nonzeros(max(switch_id_history,[],2));
                fix_qindex_history2 = switch_qindex_history(fix_id2,:);
                fix_qindex2 = max(fix_qindex_history2,[],2);  
                % append
                fix_id = [fix_id; fix_id2];
                fix_qindex = [fix_qindex; fix_qindex2];
                fix_qindex_history = [fix_qindex_history; fix_qindex_history2];
                % set flag
                fix_flag = fix_flag + 1;
                repeat = repeat*0;
            end
        end
        if fix_flag > 0
            q_opt_index(fix_id) = fix_qindex;                
        end          
        
        % update qindicator
        q_indicator = cal_indicator(param, q_opt_index); 

        % update integrals
        integ = cal_integral(const, q_indicator, eqm_cf.ExportProbQ, eqm_cf.ImportProbQ0, eqm_cf.ImportProbQ1); 
        
        % update n(q)
        nq_new = cal_nq(const, q_indicator);          
        
        % check convergence 
        err1 = max(abs(nq_new-nq));
        err2 = max(abs(eqm_cf.Pi00-Pi00));
        err3 = max(abs(eqm_cf.Pi01-Pi01));        
        err4 = max(abs(eqm_cf.Pi10-Pi10));
        err5 = max(abs(eqm_cf.Pi11-Pi11));
        err = max([err1, err2, err3, err4, err5]);
        
            % visualize progress
            scatter(const.Qgrid, nq_new, 'Marker', '.')
            xlabel('q')  
            title(strcat('iter=',num2str(iter)))
            drawnow limitrate
        
        % firm types changing quality
        switch_id = find(q_opt_index~=q_opt_index_history(:,1));
        switch_qindex = q_opt_index(switch_id);
        switch_num = size(switch_id,1);
        
        % repeat history
        for i = 1:maxHistory
            repeat(i) = repeat(i)*min(q_opt_index == q_opt_index_history(:,i)) + min(q_opt_index == q_opt_index_history(:,i));
        end                                  
        
            % print
            fprintf('iter%.0f: %.0f/%.0f types change quality', iter, switch_num, param.Nf); 
            if switch_num > 0
                if max(repeat(2:end)) > 0
                    fprintf(', repeat ='); fprintf(' %.0f', repeat); 
                end
                if switch_num <= 5
                    fprintf(', id ='); fprintf(' %.0f', switch_id);
                    fprintf(', qindex ='); fprintf(' %.0f', switch_qindex);
                end 
            elseif switch_num == 0
                fprintf(', err = %e', err);
            end
            fprintf('\n')

    end
    close
    

%% Return   
    
    eqm_cf.iter_outer      = iter;
    eqm_cf.err_outer       = err;
    if fix_flag > 0    
        eqm_cf.fix_id               = fix_id;
        eqm_cf.fix_qindex           = fix_qindex; 
        eqm_cf.fix_qindex_history   = fix_qindex_history;
    else
        eqm_cf.fix_id               = nan;
        eqm_cf.fix_qindex           = nan;         
        eqm_cf.fix_qindex_history   = nan;
    end  
    
    eqm_cf.q_index         = q_opt_index;
    eqm_cf.q_indicator     = q_indicator;
    eqm_cf.integ           = integ;  
    
    eqm_cf.nq               = nq_new;
    eqm_cf.first_nonzero    = find(eqm_cf.nq, 1, 'first');
    eqm_cf.last_nonzero     = find(eqm_cf.nq, 1, 'last');  
    eqm_cf.first_nonzero_q  = const.Qgrid(eqm_cf.first_nonzero);
    eqm_cf.last_nonzero_q   = const.Qgrid(eqm_cf.last_nonzero);  
    
    eqm_cf.density_all     = const.density_omega;  
    eqm_cf.density_exp     = const.density_omega.*full(sum(eqm_cf.ExportProbQ.*eqm_cf.q_indicator,2));
    eqm_cf.density_nonexp  = const.density_omega - eqm_cf.density_exp;

    eqm_cf.density_exp_ctn = min(eqm.density_exp, eqm_cf.density_exp);
    eqm_cf.density_nonexp_ctn = min(eqm.density_nonexp, eqm_cf.density_nonexp);
    eqm_cf.density_switcher = eqm_cf.density_all - eqm_cf.density_exp_ctn - eqm_cf.density_nonexp_ctn;    
   
    
%%
    % import probability of exporters
    eqm_cf.prob_exp_imp        = full(sum(eqm_cf.ImportProbQ1.*eqm_cf.q_indicator, 2));
    eqm_cf.prob_exp_nonimp     = full(sum((1-eqm_cf.ImportProbQ1).*eqm_cf.q_indicator, 2));
    
    % import probability of non-exporters
    eqm_cf.prob_nonexp_imp     = full(sum(eqm_cf.ImportProbQ0.*eqm_cf.q_indicator, 2));
    eqm_cf.prob_nonexp_nonimp  = full(sum((1-eqm_cf.ImportProbQ0).*eqm_cf.q_indicator, 2));
    
        % check
        check1 = eqm_cf.prob_exp_imp + eqm_cf.prob_exp_nonimp;
        check0 = eqm_cf.prob_nonexp_imp + eqm_cf.prob_nonexp_nonimp;
 
    % density of importer/non-importer
    eqm_cf.density_imp =  eqm_cf.density_exp.*eqm_cf.prob_exp_imp + eqm_cf.density_nonexp.*eqm_cf.prob_nonexp_imp;
    eqm_cf.density_nonimp = const.density_omega - eqm_cf.density_imp;

    % continuing importer etc.
    eqm_cf.density_imp_ctn = min(eqm.density_imp, eqm_cf.density_imp);
    eqm_cf.density_nonimp_ctn = min(eqm.density_nonimp, eqm_cf.density_nonimp);
    eqm_cf.density_switcher_im = eqm_cf.density_all - eqm_cf.density_imp_ctn - eqm_cf.density_nonimp_ctn;    
    
    
    
    F = eqm_cf;
 
    
end